close all; clear;

% Space-time spectrum estimation from CHAMP EEJ/EEF data
% The spatio-temporal spectra is derived following the 
% paper by Wu, Dong L., Paul B.
% Hays, Wilbert R. Skinner, 1995: A Least Squares Method for Spectral
% Analysis of Space-Time Series. J. Atmos. Sci., 52, 35013511
% Instead of LS, I use robustfit function of MATLAB
% Manoj Nair Nov 2012

%load theta

load  /Users/manojnair/projects/tides/eef_data_2000_2010.mat;


% create synthetic data for know period and wavenumber
if exist('lonnorm','var'),
    clear lonnorm;
end;

selected_lon = eej_data_2000_2010(:,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;

prd = 48*3;
wavenum = -4;

% RANDOM PHASE
% synthetic_eej = cos(2*pi*((eej_data_2000_2010(:,1)/( prd/24)) + (lonnorm'*wavenum)) + ...
% normrnd(0, 1/50, length(eej_data_2000_2010), 1) .* pi/4);

synthetic_eej = cos(2*pi*((eej_data_2000_2010(:,1)/( prd/24)) + (lonnorm'*wavenum)) + pi/4);

% significant points
% days^-1, wavenumber

[1,1;1,-3;1,-2;1,4;1,5;2,2;2-2;2,0;3,3;3,-1;4,4];

% add noise

synthetic_eej = synthetic_eej + normrnd(0, 1/40, length(eej_data_2000_2010), 1);
    
% copy the sythetic data to the actual eej array
eej_data_2000_2010(:,10) = synthetic_eej;

% # column 1: timestamp (UT) - MATLAB
% # column 2: longitude (degrees)
% # column 3: latitude (degrees)
% # column 4: local time (hours)
% # column 10: CHAMP peak current value (A/m)

% find day of the year for seasonal average

datemat = datevec(eej_data_2000_2010(:,1)); %getting date vector form date number
doyvec = dayofyear(datemat(:,1)',datemat(:,2)',datemat(:,3)'); %get DOY

% From Hermann Luehr : Here are some suggestions.
%  -  Please use the EEJ peak current density values (determined by Stefan) directly, not the eej_data_2000_2010.
%  -  Provide separate spectra for the seasons (Mai, June, July, Aug.), (Nov. Dec. Jan. Feb.),
%     (15 June - 15 Oct). The last period is centered on DE3 maximum. 4 months are in any case preferable because of the good local time coverage of CHAMP.
%  -  Is your method able to provide phase information?
%     For the full spectra you use all CHAMP years. I would suggest to separate it into the first 5 years, 
%     Aug. 2000-2005 and Aug. 2005-2010, 

days_in_prev_months = [0 31 59 90 120 151 181 212 243 273 304 334 366];
lt_min = 7;
lt_max = 17;

% All data
date_min = min(eej_data_2000_2010(:,1));
date_max = max(eej_data_2000_2010(:,1));

% Aug 2000 - Aug 2005
% date_min = min(eej_data_2000_2010(:,1));
% date_max = datenum(2005,08,1);
% 
% Aug 2005 - Aug 2010
% date_max = max(eej_data_2000_2010(:,1));
% date_min = datenum(2005,08,1);

periods = 3:1:30;
wavenum = -10:1:10;

% normalized EEJ strength
% for normalizing the effect of conductvity, I am dividing the 
% eej data with sqrt(cos((pi/12).*(LT - 12.5))))

% eef(:,10) = eef(:,10)./sqrt(cos((pi/12).*(eef(:,4) - 12.5)));
% 
% 
% 
% % delete the complex numbers 
% 
% L = (eef(:,10) == real(eef(:,10)));
% 
% eef(~L,10) = 0;

%define colormap for plots
map = [     1.0000    1.0000    1.0000
    0.9286    1.0000    1.0000
    0.8571    1.0000    1.0000
    0.7857    1.0000    1.0000
    0.7143    1.0000    1.0000
    0.6429    1.0000    1.0000
    0.5714    1.0000    1.0000
    0.5000    1.0000    1.0000
    0.4286    1.0000    1.0000
    0.3571    1.0000    1.0000
    0.2857    1.0000    1.0000
    0.2143    1.0000    1.0000
    0.1429    1.0000    1.0000
    0.0714    1.0000    1.0000
         0    1.0000    1.0000
    0.0526    1.0000    0.9474
    0.1053    1.0000    0.8947
    0.1579    1.0000    0.8421
    0.2105    1.0000    0.7895
    0.2632    1.0000    0.7368
    0.3158    1.0000    0.6842
    0.3684    1.0000    0.6316
    0.4211    1.0000    0.5789
    0.4737    1.0000    0.5263
    0.5263    1.0000    0.4737
    0.5789    1.0000    0.4211
    0.6316    1.0000    0.3684
    0.6842    1.0000    0.3158
    0.7368    1.0000    0.2632
    0.7895    1.0000    0.2105
    0.8421    1.0000    0.1579
    0.8947    1.0000    0.1053
    0.9474    1.0000    0.0526
    1.0000    1.0000         0
    1.0000    0.9444         0
    1.0000    0.8889         0
    1.0000    0.8333         0
    1.0000    0.7778         0
    1.0000    0.7222         0
    1.0000    0.6667         0
    1.0000    0.6111         0
    1.0000    0.5556         0
    1.0000    0.5000         0
    1.0000    0.4444         0
    1.0000    0.3889         0
    1.0000    0.3333         0
    1.0000    0.2778         0
    1.0000    0.2222         0
    1.0000    0.1667         0
    1.0000    0.1111         0
    1.0000    0.0556         0
    1.0000         0         0
    0.9583         0         0
    0.9167         0         0
    0.8750         0         0
    0.8333         0         0
    0.7917         0         0
    0.7500         0         0
    0.7083         0         0
    0.6667         0         0
    0.6250         0         0
    0.5833         0         0
    0.5417         0         0
    0.5000         0         0];

%% Robust Estimate of tidal signals in the CHAMPP data

% Monthly estimates of tides
for months = 1:12,

 st = days_in_prev_months(months);
 en = days_in_prev_months(months+1);

 % select data for the present month, lt and time limits
LL = doyvec >= st & doyvec <= en & eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max & eej_data_2000_2010(:,1)' > date_min & eej_data_2000_2010(:,1)' < date_max;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eej_data_2000_2010(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;


clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for i = 1:length(periods),
    for j = 1 : length(wavenum),
    
prd = periods(i);

x_m2 = [ cos(2*pi*((eej_data_2000_2010(LL,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
         sin(2*pi*((eej_data_2000_2010(LL,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
[spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eej_data_2000_2010(LL,10));
p_st(i,j,:) = stats(i,j).p;
std_err(i,j,:) = stats(i,j).se;
%spectra_ls_eej(i,j,:) = x_m2\(eej_data_2000_2010(LL,10));
end;
end;
eval(['save /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
end;

% Estimate the tides for all days

st = 1;
en = 366;

LL = doyvec >= st & doyvec <= en & eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max & eej_data_2000_2010(:,1)' > date_min & eej_data_2000_2010(:,1)' < date_max;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eej_data_2000_2010(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;



clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for i = 1:length(periods),
    for j = 1 : length(wavenum),
    
prd = periods(i);

x_m2 = [ cos(2*pi*((eej_data_2000_2010(LL,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
         sin(2*pi*((eej_data_2000_2010(LL,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
[spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eej_data_2000_2010(LL,10));
p_st(i,j,:) = stats(i,j).p;
std_err(i,j,:) = stats(i,j).se;
%spectra_ls_eej(i,j,:) = x_m2\(eej_data_2000_2010(LL,10));
end;
end;
eval(['save /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);


% Estimate for  seasons. dayofyear returns day of the year
season_st = [dayofyear(2000,05,01) dayofyear(2000,06,15) dayofyear(2000,11,01)];
season_en = [dayofyear(2000,08,30) dayofyear(2000,10,15) dayofyear(2000,02,29)];

for seasons = 1:3,

 st = season_st(seasons);
 en = season_en(seasons);


 if en > st,
     LL = doyvec >= st & doyvec <= en & eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max & eej_data_2000_2010(:,1)' > date_min & eej_data_2000_2010(:,1)' < date_max;
 else, %seasons across year borders
     LL =  ( doyvec >= st & doyvec < 366 ) | ( doyvec > 0 & doyvec <= en ) & eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max & eej_data_2000_2010(:,1)' > date_min & eej_data_2000_2010(:,1)' < date_max; % for seasons across year border
 end;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eej_data_2000_2010(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;



clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for i = 1:length(periods),
    for j = 1 : length(wavenum),
    
prd = periods(i);

x_m2 = [ cos(2*pi*((eej_data_2000_2010(LL,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
         sin(2*pi*((eej_data_2000_2010(LL,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
[spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eej_data_2000_2010(LL,10));
p_st(i,j,:) = stats(i,j).p;
std_err(i,j,:) = stats(i,j).se;
end;
end;
eval(['save /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
end;

% Estimate planetary standing waves for all days


st = 1;
en = 366;

LL = doyvec >= st & doyvec <= en & eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max & eej_data_2000_2010(:,1)' > date_min & eej_data_2000_2010(:,1)' < date_max;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eej_data_2000_2010(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;



clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for j = 1:5, % for n= 1,2,3,4,5
    

x_m2 = [ cos(2*pi*(lonnorm'*j)) ...
         sin(2*pi*(lonnorm'*j))];  
[spectra_rob_eej(j,:), stats(j) ] = robustfit(x_m2,eej_data_2000_2010(LL,10));
p_st(j,:) = stats(j).p;
std_err(j,:) = stats(j).se;
%spectra_ls_eej(i,j,:) = x_m2\(eej_data_2000_2010(LL,10));
end;

eval(['save /Users/manojnair/projects/tides/plots/' sprintf('SPW_DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

%  Monthly estimates of planetary waves SPW

clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats


for months = 1:12,

 st = days_in_prev_months(months);
 en = days_in_prev_months(months+1);

 % select data for the present month, lt and time limits
LL = doyvec >= st & doyvec <= en & eej_data_2000_2010(:,4)' >= lt_min & eej_data_2000_2010(:,4)' <= lt_max & eej_data_2000_2010(:,1)' > date_min & eej_data_2000_2010(:,1)' < date_max;
[n_data,n_days]=hist(doyvec(LL));

% make normalized latitudes between 0 (= 0) and 360 (=1)

if exist('lonnorm'),
    clear lonnorm;
end;

selected_lon = eej_data_2000_2010(LL,2);
L = selected_lon < 0;
lonnorm(L) = 360 + selected_lon(L) ;
lonnorm(~L) =  selected_lon(~L) ;
lonnorm = lonnorm/360;


clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
for j = 1:5, %n = 1,2,3,4,5
    

x_m2 = [ cos(2*pi*(lonnorm'*j)) ...
         sin(2*pi*(lonnorm'*j))];  
[spectra_rob_eej(j,:), stats(j) ] = robustfit(x_m2,eej_data_2000_2010(LL,10));
p_st(j,:) = stats(j).p;
std_err(j,:) = stats(j).se;
%spectra_ls_eej(i,j,:) = x_m2\(eej_data_2000_2010(LL,10));
end;

eval(['save /Users/manojnair/projects/tides/plots/' sprintf('SPW_DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
end;

%% Plot the data
%monthly data
close all; clear spectra_rob_eej_all_months pvalue_all_months spectra_monthly phase_monthly pvalue_monthly;

p_test_lim = 0.05;%
%p_test_lim = 100;

days_in_prev_months = [0 31 59 90 120 151 181 212 243 273 304 334 366];


nmonth = 1;

for i = [12,1:12,1], % patch december and jan months
st = days_in_prev_months(i);
en = days_in_prev_months(i+1);

% tides
    
    eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

spectra_rob_eej_all_months(nmonth,:,:,:) = spectra_rob_eej;
pvalue_all_months(nmonth,:,:,:) = p_st;



nmonth = nmonth + 1;

end;


% smooth across the months


for i = 1:12,

    A = squeeze(spectra_rob_eej_all_months(i,:,:,2) * 0.25 + spectra_rob_eej_all_months(i+1,:,:,2) * 0.5 + spectra_rob_eej_all_months(i+2,:,:,2) * 0.25);
    B = squeeze(spectra_rob_eej_all_months(i,:,:,3) * 0.25 + spectra_rob_eej_all_months(i+1,:,:,3) * 0.5 + spectra_rob_eej_all_months(i+2,:,:,3) * 0.25);
    

%     A = squeeze( spectra_rob_eej_all_months(i+1,:,:,2) ); % no smoothing
%     B = squeeze( spectra_rob_eej_all_months(i+1,:,:,3) ); % no smoothing

spectra_rob_eej_mag = sqrt(A.^2+B.^2);
spectra_rob_eej_phase = 24*((180/pi*angle(squeeze(complex(A,B))))/360);
L = squeeze(pvalue_all_months(i+1,:,:,2)) < p_test_lim & squeeze(pvalue_all_months(i+1,:,:,3)) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = NaN; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = NaN;

spectra_monthly(i,:,:) = spectra_rob_eej_mag;
phase_monthly(i,:,:) = spectra_rob_eej_phase;
pvalue_monthly(i,:,:) = max ( squeeze(p_st(:,:,2)), squeeze(p_st(:,:,3)));
end;

% plot Non-migrating dirunal

fig=figure(1);
set(fig,'Position', [670         825        1007         265]); % set paper size
[y2,ia2,ib2] = intersect(wavenum,[-4,-3,-2,-1,0,2,3,4,5]); % select wavenumbers
[y1,ia1,ib1] = intersect(periods,24); % select period
image_data = (squeeze(spectra_monthly(:,ia1,ia2)))'; % get the data 
LLL = isnan(image_data(:)); % NaN's indicate points that are not statistically significant
LL = isnan(image_data); % column vector for text plotting, matrix for image
% print data matrix

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','wt');

fprintf(fid,'Nonmigrating dirunal Amplitude\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
                fprintf(fid,'%6.2f ', image_data(kk,kkk)); 
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

image_data(LL) = nanmin(image_data(:))-(nanmax(image_data(:))-nanmin(image_data(:)))/10; % trick to get high P points to be white
imagesc(datenum(2000,1:12,0),1:length(ia2),image_data );
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 2';' 3';' 4';' 5']);
title(['Non migrating dirunal' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
colormap(map);
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;

text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');
close all;






% plot Non-migrating semi-dirunal

fig=figure(2);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,12);
[y2,ia2,ib2] = intersect(wavenum,[-3,-2,-1,0,1,3,4,5,6]);
%imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(spectra_monthly(:,ia1,ia2)))');
image_data = (squeeze(spectra_monthly(:,ia1,ia2)))';
LLL = isnan(image_data(:));
LL = isnan(image_data);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Non-migrating semi-dirunal Amplitude\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
                fprintf(fid,'%6.2f ', image_data(kk,kkk)); 
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

image_data(LL) = nanmin(image_data(:))-(nanmax(image_data(:))-nanmin(image_data(:)))/10;
imagesc(datenum(2000,1:12,0),1:length(ia2),image_data );
set(gca,'YTick',1:length(ia2));
%set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 1';' 3';' 4']);
set(gca,'YTickLabel',['-3';'-2';'-1';' 0';' 1';' 3';' 4';' 5';' 6']);
title(['Non migrating semi-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
colormap(map);
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');

close all;

% plot Non-migrating teri-dirunal
fig=figure(3);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,8);
[y2,ia2,ib2] = intersect(wavenum,[-2,-1,0,1,2,4,5,6,7]);
%imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(spectra_monthly(:,ia1,ia2)))');
image_data = (squeeze(spectra_monthly(:,ia1,ia2)))';
LLL = isnan(image_data(:));
LL = isnan(image_data);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Non-migrating teri-dirunal Amplitude\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
                fprintf(fid,'%6.2f ', image_data(kk,kkk)); 
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

image_data(LL) = nanmin(image_data(:))-(nanmax(image_data(:))-nanmin(image_data(:)))/10;
imagesc(datenum(2000,1:12,0),1:length(ia2),image_data );
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-2';'-1';' 0';' 1';' 2';' 4';' 5';' 6';' 7']);
title(['Non migrating teri-dirunal' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
colormap(map);
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');


saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');
close all;

% plot Non-migrating q-dirunal
fig=figure(4);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,6);
[y2,ia2,ib2] = intersect(wavenum,[-1,0,1,2,3,5,6,7,8]);
image_data = (squeeze(spectra_monthly(:,ia1,ia2)))';
LLL = isnan(image_data(:));
LL = isnan(image_data);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Non-migrating q-dirunal Amplitude\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
                fprintf(fid,'%6.2f ', image_data(kk,kkk)); 
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

image_data(LL) = nanmin(image_data(:))-(nanmax(image_data(:))-nanmin(image_data(:)))/10;
imagesc(datenum(2000,1:12,0),1:length(ia2),image_data );
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-1';' 0';' 1';' 2';' 3';' 5';' 6';' 7';' 8']);
title(['Non migrating q-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
colormap(map);
h = colorbar; ylabel(h, 'EEJ amplitude mA/m');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');


saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');



% plot phase of the non-migrating tides

close all;
% plot Non-migrating dirunal
fig=figure(1);
set(fig,'Position', [670         825        1007         265]);
[y2,ia2,ib2] = intersect(wavenum,[-4,-3,-2,-1,0,2,3,4,5]);
image_data = squeeze(phase_monthly(:,22,ia2))';
LLL = isnan(image_data(:));
LL = isnan(image_data);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Non-migrating dirunal phase\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
                fprintf(fid,'%6.2f ', image_data(kk,kkk)); 
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

image_data(LL) = -12;
imagesc(datenum(2000,1:12,0),1:length(ia2),image_data );
%imagesc(datenum(2000,1:12,15),1:length(ia2), squeeze(phase_monthly(:,22,ia2))');
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 2';' 3';' 4';' 5']);
title(['Non migrating dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');

h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');

saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Dirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;
% plot Non-migrating semi-dirunal
fig=figure(2);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,12);
[y2,ia2,ib2] = intersect(wavenum,[-3,-2,-1,0,1,3,4,5,6]);
%imagesc(datenum(2000,1:12,15),1:length(ia2), (squeeze(phase_monthly(:,ia1,ia2)))');
image_data = (squeeze(phase_monthly(:,ia1,ia2)))';

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Non-migrating semi-dirunal phase\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
        fprintf(fid,'%6.2f ', image_data(kk,kkk));
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

LLL = isnan(image_data(:));
LL = isnan(image_data);
image_data(LL) = -12;
imagesc(datenum(2000,1:12,0),1:length(ia2),image_data );
set(gca,'YTick',1:length(ia2));
%set(gca,'YTickLabel',['-4';'-3';'-2';'-1';' 0';' 1';' 3';' 4']);
set(gca,'YTickLabel',['-3';'-2';'-1';' 0';' 1';' 3';' 4';' 5';' 6']);
title(['Non migrating semi-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
%colormap(map);
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');


saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Semidirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;
% plot Non-migrating teri-dirunal
fig=figure(3);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,8);

[y2,ia2,ib2] = intersect(wavenum,[-2,-1,0,1,2,4,5,6,7]);
image_data = (squeeze(phase_monthly(:,ia1,ia2)))';
LLL = isnan(image_data(:));
LL = isnan(image_data);


fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Non-migrating  teri-dirunal phase\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
        fprintf(fid,'%6.2f ', image_data(kk,kkk));
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

image_data(LL) = -12;
imagesc(datenum(2000,1:12,0),1:length(ia2),image_data );
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-2';'-1';' 0';' 1';' 2';' 4';' 5';' 6';' 7']);
title(['Non migrating teri-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
%colormap(map);
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');



saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Teridirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;
% plot Non-migrating q-dirunal
fig=figure(4);
set(fig,'Position', [670         825        1007         265]);
[y1,ia1,ib1] = intersect(periods,6);
[y2,ia2,ib2] = intersect(wavenum,[-1,0,1,2,3,5,6,7,8]);
image_data = (squeeze(phase_monthly(:,ia1,ia2)))';
LLL = isnan(image_data(:));
LL = isnan(image_data);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Non-migrating q-dirunal phase\n');
for kk = 1:size(image_data,1),
    for kkk = 1:size(image_data,2),
        fprintf(fid,'%6.2f ', image_data(kk,kkk));
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

image_data(LL) = -12;
imagesc(datenum(2000,1:12,0),1:length(ia2), image_data);
set(gca,'YTick',1:length(ia2));
set(gca,'YTickLabel',['-1';' 0';' 1';' 2';' 3';' 5';' 6';' 7';' 8']);
title(['Non migrating Q-dirunal ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))]);
ylabel('wavenumber');
%colormap(map);
h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:length(ia2));
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(round(image_data(:)*100)/100)));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');

set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Qdirunal_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');


close all;


% Plot migrating tides using spectra_monthly and phase_monthly arrays from the above cell


[y1,ia1,ib1] = intersect(periods,[3,4,6,8,12,24]);
[y2,ia2,ib2] = intersect(wavenum,[1,2,3,4,6,8]);


j = 1;
for i = length(ia1):-1:1,
    
mig_mat(j,:) = (squeeze(spectra_monthly(:,ia1(i),ia2(j))));
mig_mat_phase(j,:) = (squeeze(phase_monthly(:,ia1(i),ia2(j))));
j = j+1;
end;

fig=figure(1);
set(fig,'Position', [670         825        1007         265]);

LLL = isnan(mig_mat(:));
LL = isnan(mig_mat);
mig_mat(LL) = nanmin(mig_mat(:))-(nanmax(mig_mat(:))-nanmin(mig_mat(:)))/100;

imagesc(datenum(2000,1:12,0),1:6, log10(mig_mat));

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Migrating Amplitude\n');

for kk = 1:size(mig_mat,1),
    for kkk = 1:size(image_data,2),
        fprintf(fid,'%6.2f ', log10(mig_mat(kk,kkk)));
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

colormap(map);
set(gca,'YTick',1:6);
set(gca,'YTickLabel',['1,1';'2,2';'3,3';'4,4';'6,6';'8,8']);
title(['Amplitude of migrating waves ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))] );
ylabel('n,s');
h = colorbar; ylabel(h, 'EEJ amplitude log (mA/m) ');
datetick('x','mmm','keeplimits');

% plot numbers in the cell

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:6);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(log10(mig_mat(:)))));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');


set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');


close all;
% phase of the tides
fig=figure(2);
set(fig,'Position', [670         825        1007         265]);
LLL = isnan(mig_mat_phase(:));
LL = isnan(mig_mat_phase);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Migrating Phase\n');
for kk = 1:size(mig_mat,1),
    for kkk = 1:size(mig_mat_phase,2),
        fprintf(fid,'%6.2f ', mig_mat_phase(kk,kkk));
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

mig_mat_phase(LL) = -13;
imagesc(datenum(2000,1:12,0),1:6, (mig_mat_phase));
set(gca,'YTick',1:6);
set(gca,'YTickLabel',['1,1';'2,2';'3,3';'4,4';'6,6';'8,8']);
title(['Phase of migrating waves ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))] );
ylabel('n,s');

% plot numbers in the cell

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:6);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(mig_mat_phase(:))));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');

h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Migrating_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;


% Plot the seasons data
season_st = [dayofyear(2000,05,01) dayofyear(2000,06,15) dayofyear(2000,11,01)];
season_en = [dayofyear(2000,08,30) dayofyear(2000,10,15) dayofyear(2000,02,29)];
fig = figure(1);
for seasons = 1:3,

 st = season_st(seasons);
 en = season_en(seasons);


 eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2);
spectra_rob_eej_phase = 180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3))));
L = p_st(:,:,2) < p_test_lim & p_st(:,:,3) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;
subplot(length(season_st),2,(seasons-1)*2 + 1);
imagesc( wavenum, periods,log10(spectra_rob_eej_mag));
colormap(map);
set(gca,'FontSize',10);
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));
axis([-6 10 3 26]);
caxis([0,2]);
h = colorbar; ylabel(h, 'EEJ amplitude log(mA/m)');;
subplot(length(season_st),2,(seasons-1)*2 + 2);
imagesc( wavenum, periods,24 * ( spectra_rob_eej_phase./360 ));
%axis([-6 10 3 26]);
axis([-6 10 2.5 26.5]);
caxis([-12,12]);
set(gca,'FontSize',10);
 h = colorbar; ylabel(h, 'Phase hours');
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));

end;
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Seasons_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/Seasons_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'png');

close all;
% plot ALL data

st = 1;
en = 366;

 eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2);
spectra_rob_eej_phase = 180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3))));
L = p_st(:,:,2) < p_test_lim & p_st(:,:,3) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = NaN; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;
figure(1);
imagesc( wavenum, periods,log10(spectra_rob_eej_mag));
colormap(map);
set(gca,'FontSize',16);
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));
%axis([-6 10 2 26]);
axis([-6 10 2.5 26.5]);
caxis([0,2]);
h = colorbar; ylabel(h, 'EEJ amplitude log(mA/m)');
xlabel('wavenumber');
ylabel('period (hours)');

saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Amplitude_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Amplitude_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'png');

figure(2);
LLL = isnan(spectra_rob_eej_phase);
spectra_rob_eej_phase(LLL) = -180;
imagesc( wavenum, periods,24 * ( spectra_rob_eej_phase./360 ));
%axis([-6 10 3 26]);
axis([-6 10 2.5 26.5]);
caxis([-12,12]);
set(gca,'FontSize',16);
h = colorbar; ylabel(h, 'Phase in hours');
title( sprintf('DOY %d - %d,LT %d - %d,p-limit %4.2f,year %s-%s',st,en, lt_min,lt_max, p_test_lim,datestr(date_min,10),datestr(date_max,10)));
set(fig, 'PaperPositionMode', 'manual');
xlabel('wavenumber');
ylabel('period (hours)');
colormap(colormap_map_phase);
colormap_map_phase = colormap;
colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Phase_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/ALL_Phase_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10))],'png');

close all;


%% Plot Standing Planetary Waves Amplitude Monthly

clear spectra_rob_eej_all_months pvalue_all_months spectra_monthly phase_monthly pvalue_monthly;

nmonth = 1;

for i = [12,1:12,1], % patch december and jan months
st = days_in_prev_months(i);
en = days_in_prev_months(i+1);

% tides
    
    eval(['load /Users/manojnair/projects/tides/plots/' sprintf('SPW_DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);

spectra_rob_eej_all_months(nmonth,:,:) = spectra_rob_eej;
pvalue_all_months(nmonth,:,:) = p_st;



nmonth = nmonth + 1;

end;


% smooth across the months


for i = 1:12,

    A = squeeze(spectra_rob_eej_all_months(i,:,2) * 0.25 + spectra_rob_eej_all_months(i+1,:,2) * 0.5 + spectra_rob_eej_all_months(i+2,:,2) * 0.25);
    B = squeeze(spectra_rob_eej_all_months(i,:,3) * 0.25 + spectra_rob_eej_all_months(i+1,:,3) * 0.5 + spectra_rob_eej_all_months(i+2,:,3) * 0.25);
    

%     A = squeeze( spectra_rob_eej_all_months(i+1,:,:,2) ); % no smoothing
%     B = squeeze( spectra_rob_eej_all_months(i+1,:,:,3) ); % no smoothing

spectra_rob_eej_mag = sqrt(A.^2+B.^2);
spectra_rob_eej_phase = 24*((180/pi*angle(squeeze(complex(A,B))))/360);
L = squeeze(pvalue_all_months(i+1,:,2)) < p_test_lim & squeeze(pvalue_all_months(i+1,:,3)) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = NaN; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = NaN;

spectra_monthly(i,:) = spectra_rob_eej_mag;
phase_monthly(i,:) = spectra_rob_eej_phase;
pvalue_monthly(i,:) = max ( squeeze(p_st(:,2)), squeeze(p_st(:,3)));
end;

mig_mat = spectra_monthly';
mig_mat_phase = phase_monthly';

fig=figure(1);
set(fig,'Position', [670         825        1007         265]);

LLL = isnan(mig_mat(:));
LL = isnan(mig_mat);
mig_mat(LL) = nanmin(mig_mat(:))-(nanmax(mig_mat(:))-nanmin(mig_mat(:)))/100;

imagesc(datenum(2000,1:12,0),1:5, mig_mat);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Planetary Standing Wave Amplitude\n');

for kk = 1:size(mig_mat,1),
    for kkk = 1:size(mig_mat,2),
        fprintf(fid,'%6.2f ', mig_mat(kk,kkk));
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

colormap(map);
set(gca,'YTick',1:5);
set(gca,'YTickLabel',['1';'2';'3';'4';'5']);
title(['Amplitude of planetary standing waves ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))] );
ylabel('n');
h = colorbar; ylabel(h, 'EEJ amplitude (mA/m) ');
datetick('x','mmm','keeplimits');

% plot numbers in the cell

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:5);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(mig_mat(:))));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');


set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/SPW_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/SPW_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Amplitude'],'png');


%close all;


% phase of the planetary standing waves

fig=figure(2);

set(fig,'Position', [670         825        1007         265]);
LLL = isnan(mig_mat_phase(:));
LL = isnan(mig_mat_phase);

fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
fprintf(fid,'Planetary standing waves Phase\n');
for kk = 1:size(mig_mat_phase,1),
    for kkk = 1:size(mig_mat_phase,2),
        fprintf(fid,'%6.2f ', mig_mat_phase(kk,kkk));
    end;
    fprintf(fid,'\n');
end;
fclose(fid);

%mig_mat_phase(LL) = -14;
imagesc(datenum(2000,1:12,0),1:5, (mig_mat_phase));
set(gca,'YTick',1:5);
set(gca,'YTickLabel',['1';'2';'3';'4';'5']);
title(['Phase of Planetary Standing Waves ' sprintf(' LT %d-%d p-limit %4.2f year %s-%s',lt_min,lt_max,p_test_lim,datestr(date_min,10),datestr(date_max,10))] );
ylabel('n');

% plot numbers in the cell

[xlbl, ylbl] = meshgrid(datenum(2000,1:12,0),1:5);
%# create cell arrays of number labels
lbl = strtrim(cellstr(num2str(mig_mat_phase(:))));
for ii = 1:length(lbl), index_i = find(lbl{ii} =='.');  lbl{ii}(index_i+3:end) = []; end;
colormap_map_phase = colormap;
%colormap_map_phase(1,:) = [1,1,1];
colormap(colormap_map_phase);
text(xlbl(~LLL), ylbl(~LLL), lbl(~LLL),'color','k',...
    'HorizontalAlignment','center','VerticalAlignment','middle');

h = colorbar; ylabel(h, 'Phase hours');
datetick('x','mmm','keeplimits');
set(fig, 'PaperPositionMode', 'manual');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/SPW_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'fig');
saveas(gcf, ['/Users/manojnair/projects/tides/plots/SPW_Month_' sprintf('LT_%d_%d_Year_%s_%s',lt_min,lt_max,datestr(date_min,10),datestr(date_max,10)) 'Phase'],'png');

close all;


% End standing planetary waves

%% Print the results as a table
p_test_lim = 100; % Just to make sure that we don't print garbage !
fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');

fprintf(fid,'Period (Hour),Wavenumber  A(mA/m) B(mA/m) Std_Err_A Std_Err_B Amp(mA/m) Phase(Hours)\n');
season_st = [dayofyear(2000,05,01) dayofyear(2000,06,15) dayofyear(2000,11,01)];
season_en = [dayofyear(2000,08,30) dayofyear(2000,10,15) dayofyear(2000,02,29)];

for seasons = 1:3,

 st = season_st(seasons);
 en = season_en(seasons);
 
fprintf(fid,'DOY_%d_%d_LT_%d_%d_T_%s_%s\n',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29));
eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2);
spectra_rob_eej_phase = 24*((180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3)))))/360);
L = p_st(:,:,2) < p_test_lim & p_st(:,:,3) < p_test_lim; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;

[y1,ia1,ib1] = intersect(periods,[24,12,8,6,4]);
[y2,ia2,ib2] = intersect(wavenum,[-7:7]);

for i = 1:length(ia1),
    for j = 1:length(ia2),
        
        if p_st(ia1(i),ia2(j),2) <0.05 & p_st(ia1(i),ia2(j),3) <0.05,
            fprintf(fid,'%2d,%2d   %6.2f %6.2f %6.2f %6.2f %6.2f %5.2f\n', ...
                periods(ia1(i)),wavenum(ia2(j)), ...
                spectra_rob_eej(ia1(i),ia2(j),2), ...
                spectra_rob_eej(ia1(i),ia2(j),3), ...
                std_err(ia1(i),ia2(j),2), ...
                std_err(ia1(i),ia2(j),3), ...
                spectra_rob_eej_mag(ia1(i),ia2(j)), ...
                spectra_rob_eej_phase(ia1(i),ia2(j)));
        end;
    end;
end;
                
    
end;


fclose(fid);


% print the results for all data
p_test_lim = 0.05; % Just to make sure that we don't print garbage !
fid = fopen('/Users/manojnair/projects/tides/plots/ionospheric_results.txt','at');
st = 1;
en = 366;

 eval(['load /Users/manojnair/projects/tides/plots/' sprintf('DOY_%d_%d_LT_%d_%d_T_%s_%s',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29)) '_ROB_SPECTRA' ...
    ' spectra_rob_eej p_st std_err st en periods wavenum n_data n_days']);
fprintf(fid,'DOY_%d_%d_LT_%d_%d_T_%s_%s\n',st,en,lt_min,lt_max,datestr(date_min,29),datestr(date_max,29));
spectra_rob_eej_mag = sqrt(spectra_rob_eej(:,:,2).^2+spectra_rob_eej(:,:,3).^2);
spectra_rob_eej_phase = 24*((180/pi*angle(squeeze(complex(spectra_rob_eej(:,:,2),spectra_rob_eej(:,:,3)))))/360);
L = p_st(:,:,2) < 0.05 & p_st(:,:,3) < 0.05; % Get the estimates for which both the real and imagianry parts are statistically significat
spectra_rob_eej_phase(~L) = 0; % Assign 0 to insigificant grid points
spectra_rob_eej_mag(~L) = 0;

[y1,ia1,ib1] = intersect(periods,[24,12,8,6,4]);
[y2,ia2,ib2] = intersect(wavenum,[-7:7]);

for i = 1:length(ia1),
    for j = 1:length(ia2),
        
        if p_st(ia1(i),ia2(j),2) <0.05 & p_st(ia1(i),ia2(j),3) <0.05,
            fprintf(fid,'%2d,%2d   %6.2f %6.2f %6.2f %6.2f %6.2f %5.2f\n', ...
                periods(ia1(i)),wavenum(ia2(j)), ...
                spectra_rob_eej(ia1(i),ia2(j),2), ...
                spectra_rob_eej(ia1(i),ia2(j),3), ...
                std_err(ia1(i),ia2(j),2), ...
                std_err(ia1(i),ia2(j),3), ...
                spectra_rob_eej_mag(ia1(i),ia2(j)), ...
                spectra_rob_eej_phase(ia1(i),ia2(j)));
        end;
    end;
end;
fclose(fid);


% 
% %%
% %testing whether ordered longitude is necessary
% 
% [y,ia] = sort(eej_data_2000_2010(:,2));
% 
% x_lon = [ cos(2*pi*((lonnorm(ia)'*2))) ...
%           sin(2*pi*((lonnorm(ia)'*2)))];  
% 
% spectra_ln = x_lon\eej_data_2000_2010(ia,6);
% 
% % now invert without ordering
% 
% x_lon_n = [ cos(2*pi*((lonnorm'*3.5))) ...
%           sin(2*pi*((lonnorm'*3.5)))];  
% 
% spectra_ln_n = x_lon_n\eej_data_2000_2010(:,6);
% 
% %inference. There is no difference between ordered or non-ordered inversion
% 
% % finding out the wave structure of M2  tides
% 
% 
% periods = 12.415:0.001:12.43 ; 
% wavenum = -8:1:8;
% clear spectra_ls spectra_rob_eej spectra_ls_eej
% for i = 1:length(periods),
%     for j = 1 : length(wavenum),
%     
% prd = periods(i);
% 
% x_m2 = [ cos(2*pi*((eej_data_2000_2010(:,1)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
%          sin(2*pi*((eej_data_2000_2010(:,1)/( prd/24)) + (lonnorm'*wavenum(j))))];  
% [spectra_rob_eej(i,j,:), stats(i,j) ] = robustfit(x_m2,eej_data_2000_2010(:,10));
% %spectra_ls_eej(i,j,:) = x_m2\(eej_data_2000_2010(LL,10));
% end;
% end;
% 

%% Joint CHAMP eej_data_2000_2010 and Jicamarca ISR eej_data_2000_2010 analysis
% The idea is to include the Icamarca ISR data to the inversion. This
% should remove the "nigh-time" ambiguity in the wave fitting. Inference.
% By including the ISR data, the wavemodes are no longer discrete, but
% continous signals, with the highest at the migrating tides mode. For
% example, the relative differences between DE3 and SE2 is not very clear
% now. However, one should also look whether Jicamarca data is overwhelming the
% inversion process. For example, when I limited jicamarca data to just
% 5000 samples, the spectral image became much similar to the one I got
% without including Jicamarca data.

% load /data/backup/mnair/ace/jcamarca_isr_fejer.mat eej_data_2000_2010;
% jeef = eef;
% load /data/backup/mnair/longp/eef_data_2000_2009.mat eef;
% 
% 
% jeef(:,1) = jeef(:,1) + datenum(2000,1,1);
% 
% % running mean of Jicamarca ISR 
% 
% jic_eef = runmean(jeef(1:5000,6),1); % 
% 
% 
% % add jicamarca data to champ data
% 
% t = [eef(:,1)' jeef(1:5000,1)']';
% y = [eef(:,6)' jic_eef']';
% l = [eef(:,2)' ones([1,length(jic_eef)])*-76.8];
% lt = [eef(:,4)' jeef(1:5000,4)'];
% 
% lt_min = 7;
% lt_max = 17;
% date_min = datenum(2000,08,1);
% date_max = datenum(2005,07,31);
% p_test_lim = 0.05;
% periods = 3:1:30;
% wavenum = -10:1:10;
% 
% 
% LL =  t > date_min & t < date_max;
% 
% [n_data,n_days]=hist(doyvec(LL));
% 
% % make normalized latitudes between 0 (= 0) and 360 (=1)
% 
% if exist('lonnorm'),
%     clear lonnorm;
% end;
% 
% selected_lon = l(LL);
% L = selected_lon < 0;
% lonnorm(L) = 360 + selected_lon(L) ;
% lonnorm(~L) =  selected_lon(~L) ;
% lonnorm = lonnorm/360;
% 
% 
% clear spectra_ls spectra_rob_eej spectra_ls_eej p_st std_err stats
% for i = 1:length(periods),
%     for j = 1 : length(wavenum),
%     
% prd = periods(i);
% 
% x_m2 = [ cos(2*pi*((t(LL)/( prd/24)) + (lonnorm'*wavenum(j)))) ...
%          sin(2*pi*((t(LL)/( prd/24)) + (lonnorm'*wavenum(j))))];  
% %[spectra_rob_eej_jic(i,j,:), stats(i,j) ] = robustfit(x_m2,y(LL));
% % p_st(i,j,:) = stats(i,j).p;
% % std_err(i,j,:) = stats(i,j).se;
% spectra_ls_eej_jic2(i,j,:) = x_m2\(y(LL));
% end;
% end;
% 

